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We extend the work of A. Ciaffaglione and P. di Gianantonio on mechanical verification 
of algorithms for exact computation on real numbers, using infinite streams of digits 
implemented as a co-inductive type. Four aspects are studied: the first aspect concerns 
the proof that digit streams can be related to axiomatized real numbers when they are 
already present in the proof system. The second aspect re-visits the definition of an 
addition function, looking at techniques to let the proof search engine perform the 
effective construction of an algorithm that is correct by construction. The third aspect 
concerns the definition of a function to compute affine formulas with positive rational 
coefficients. This is an example where we need to combine co-recursion and recursion. 
The fourth aspect concerns the definition of a function to compute series, with an 
application on the series that is used to compute Euler's number e. All these 
experiments should be reproducible in any proof system that supports co-inductive 
types, co-recursion and general forms of terminating recursion. We used the COQ system 
Powek et al., 1993| [Bertot and Casteran, 2004||Gimenez, 1994^ . 

1. Introduction 

Several proof systems provide data-types to describe real numbers, together with basic op- 
erations and theorems giving an ordered, complete, and archimedian field ( [Harrison, 1996| 
[Harrison, 1998| [Mayero, 2001| ). In the COQ system, several approaches have been taken; 
depending on whether developers wanted to adhere to pure constructive mathematics or 
more classical approaches. In the classical approach, the type of real numbers is merely 
"axiomatized" , the existence of the type and the elementary operations is assumed and 
the properties of these operations are asserted as axioms. This approach has been used 
extensively to provide a large collection of results, going all the way to the description of 
trigonometric functions, calculus and the like. However, because the type of real numbers 
is axiomatized, there is no "physical representation of numbers" and the basic operations 
correspond to no algorithm. 

In an alternative approach, a type of constructive numbers may be defined as a data- 
type and the basic operations may be described as algorithms manipulating elements of 
this data-type. This approach is used for instance in C-CoRN ( [Cruz-Filipe et al., 2004[ ). 
A. Ciaffaglione and P. di Gianantonio ( [Ciaffaglione and di Gianantonio, 2000[ ) showed 
that a well-known representation of real numbers as infinite sequences of redundant digits 
could easily be implemented inside theorem provers with co-inductive types. We say the 
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digits are redundant because several representation are possible for every number. In 
the case of | |Ciaffaglione and di Gianantonio, 2000| ) the representation is simply inspired 
by the usual binary representation of fractional numbers and made more redundant by 
adding the possibility to use a negative digit. A. Ciaffaglione and P. di Gianantonio then 
provide addition, multiplication, and show that these operations enjoy the properties 
that are expected from a constructive field of real numbers. 

In our approach there is also an extra digit, but its meaning is intermediary between the 
existing and 1. Edalat and Potts ( |Edalat and Potts, 1998| ) show that both approaches 
are special cases of a more general family of representations based on linear fractional 
transformations. 

Once we have chosen the way to represent real numbers as a data-type, we proceed by 
establishing a relation between this data-type and the axiomatized type of real numbers. 
This is a departure from the prescription of pure constructive mathematics, because we 
rely on the non-constructive axioms of that theory to state the correctness of our algo- 
rithms. Still, we also describe the relation between our representation and constructive 
approaches to real numbers, by showing how infinite sequences of digits can be produced 
from constructive views of Dedekind cuts or Cauchy sequences. 

Once we have provided the basic data-type and its relation with the axiomatized theory 
of real numbers, we proceed by defining an addition function. We rely on proof search 
tools to construct the addition algorithm: we only provide guidelines for the construction 
of the algorithm, without actually describing all 25 cases in the function. The correctness 
proof then consists in showing that there is a morphism between the data-type and the 
axiomatized type. Our contribution in this part is to show how to use the proof search 
engine to construct a well formed addition function. 

We then focus on afRne formulas combining two real values with rational coefficients. 
For these more general operations, we need to combine co-recursion and well-founded 
recursion. We show that the function responsible for producing the infinite stream of 
digits representing the result can be decomposed in two recursive functions. One of the 
functions is a guarded co- recursive function as proposed in ( |Gimenez, 1994| |, the second 
function is a well-founded recursive function as in ( [Nordstrom, 1988| ). Each function 
satisfies a different form of constraint: the co-recursive function does not need to be 
terminating, but it must produce at least a digit at each recursive call, while the well- 
founded function does not need to produce a digit at each recursive call, but it must 
terminate. Our main contribution in this part is to show that some functions that appear 
at first sight to be outside the expressive power of guarded co-recursion can actually be 
modeled and proved correct. 

In a fourth part, we focus on constructively converging infinite sums. We show how 
to avoid having to consider the infinity of terms that are parts of the sum. We exhibit a 
framework that can be re-used from one series to the other. As applications, we show how 
to compute the infinite stream representing Euler's number e and to multiply two real 
numbers represented as infinite streams of digits. In particular, the algorithm we obtain 
can be executed directly using the reduction mechanism provided in COQ to compute e 
to a great precision in a reasonable time. 
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Our account stops here, although the experiments described in this paper seem to open 
the door for a more complete study of real functions, especially analytic functions. 

2. Related work 

For numerical computation, real numbers are usually represented as approximates using 
floating point numbers. These floating point numbers are composed of a mantissa and 
an exponent, so that the value of the least significant bit in the mantissa varies with 
the exponent. Still both the exponent and the mantissa have a fixed size, so that there 
is only a finite number of floating point numbers and real values must be rounded to 
find the closest floating point numbers. Floating-point based computations are thus only 
approximative and errors stemming from successive rounding operations may accumulate 
to the point that some computations can become grossly wrong ( IMuller, 1997| ). 

In spite of their limitations, Floating-point numbers are used extensively: most proces- 
sors directly provide an implementation of the elementary operations (addition, subtrac- 
tion, multiplication, division) according to a standard that gives a precise mathematical 
meaning to the rounding operations ( |IEEE, 1987| ). This standard provides the basis to 
implement computations with a guaranteed precision ( |Muller, 1997||nanrot et al., 2000| ), 
sometimes with correctness proofs that can be verified with the help of computer- 
aided proof tools ( |Daumas et al., 200 1| [Harrison, 2000| |Russinoff, 1999i |Boldo, 2004| ). 
In particular, some approaches, named expansions make it possible to increase dras- 
tically the number of representable numbers by extending the length of the mantissa 
< |Boldo et al., 2002| ). 

An alternative to floating point and rounding concentrates on data-structures that sup- 
port exact computation. Among the possible approaches, the best well-known are based 
on continued fractions fCospe r, 1972||Vuillemin, 1990| | or representations with mantissas 
that grow as needed ( Menissier-Morain, 1995| ). In the latter case, the representation is 
very close to the floating-point representations with rounding modes. One way to under- 
stand this "growing mantissa" is to view it as an infinite sequence of digits, where only a 
finite prefix is known at any time. An introduction to exact real arithmetics can be found 
in ( |Edalat and Potts, 1998| ) . Implementations are provided in the setting of conventional 
programming languages ( |Lambov, 2005| [Miiller, 2005| ), or in the setting of functional 
programming HBoehm et al., 1986| [Menissier-Morain, 1994| [Bauer et al., 2002D . 

Formal proofs about computations on infinite data-structures are a privileged ground 
for the use of co-inductive types ( [Coquand, 1993||Gimenez, 1994| ). First experiments on 
the topic of exact real number computations using co-inductive types were performed 
by Ciaffaglione and di Gianantonio jCiaffaglione and di Gianantonio, 2000| ) who showed 
that one could represent infinite sequences of digits with co-inductive types and the basic 
operations of arithmetics (addition, multiplication, comparison) with simple co-recursive 
functions, as long as the set of digits was extended to allow for enough redundancy. 
Niqui (Niqui, 2004| ) also studies the problems of modeling real number arithmetic for 
use in formal proofs, providing a single point of view to account for both continued 
fractions and infinite sequences of digits. Our approach is very similar to Ciaffaglione 
and di Gianantonio's, it only differs in the collection of digits that we consider. 
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The example of combination of co-recursive functions and well-founded recursive func- 
tions that we exhibit in our treatment of affine formulas is related to work by di Gianan- 
tonio and Miculan ( |di Gianantonio and Miculan, 2003| ) and our own work on partial 
CO- recursive functions ( |Bertot, 2005| ). 

Concerning the computation of series, we are aware through personal communica- 
tion that computations of e had already been done using the real numbers as they are 
formalized in the C-CoRN library ( |Cruz-Filipe et al., 2004| ); it seems that co- inductive 
presentations make it possible to achieve a higher efficiency. While the C-CoRN library 
aims at providing a comprehensive study of constructive mathematics, our work aban- 
dons some of the foundations of constructive mathematics: we let our logical reasoning 
depend on non-constructive arguments, although we do provide algorithms working on 
concrete data-structure to represent real numbers and the basic operations. We believe 
the algorithms we develop will be useful even in the context of constructive mathematics, 
but the proofs of correctness will probably have to be re-done. 



3. Redundant digit representation for real numbers 

We are all used to the notation with a decimal point to represent real numbers. For 
instance, we usually write a number between and 1 as a string of the form 0.1354647 • • • 
and we know that the sequence of digits must be infinite for some numbers, actually all 
those that are not of the form where a and b are positive integers. It is a bit less 
natural, but still easy to understand, that all numbers can be represented by an infinite 
sequence: for those that have a finite representation, it suffices to add an infinite sequence 
of zeros. Moreover the number 1 can also be represented by the sequence 0.999 ■ • ■ 

When we know a prefix of one of these infinite sequences, we actually know the number 
that is represented up to a certain precision. If the prefix has length p, we actually 
know precisely the bounds of an interval of width that contains the number. We are 
accustomed to reasoning with these prefixes of infinite sequences and we expect tools to 
return correct prefixes of an operation's output when this operation has been fed with 
correct prefixes for the inputs. 

In the conventional representation, the number ten plays a special role: it is the base. 
We can change the base and use digits that are between and the base. For instance, 
we can use two as the base, so that the digits are only and 1. The number ^ can then 
be represented by the sequence 0.1000 ■ • ■ and the number 1 can be represented by the 
sequence 0.1111 • • •. For a sequence 0.did2 ■ ■ •, the number being represented is: 



E 

i=l 



the following equalities hold: 



0.0. = ^ 
2 



0.1s 



O.s + 1 



2 

A prefix with p digits gives an interval of width ^ that contains the number represented 
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by the infinite sequence. In the rest of this paper, we wiU carry on using base 2 (but the 
set of digits will change). 

For the computation of basic operations, the base-2 conventional representation is not 
really well adapted. Here is an example that exhibits the main flaw of this representation. 
The numbers ^ and g add up to give i. However, the infinite sequences for these numbers 
are given in the following equations: 

0.01010101- •• 
0.00101010••• 

0.10000000 • • • = 0.01111111 • • • 

The following reasoning steps justify the first equation: 

oo 

0.01010101 . • . = y — = ^ - 1 = - 

1=1 4 

Similar justifications can be used for the other equations. If w is a prefix of 0.0101 ■ • •, 
then w is also the prefix of all numbers between wOO ■ ■ ■ and will • • •. These two numbers 
are rational numbers of the form ^ and neither can be equal to ^. Thus, we actually 
have an interval of possible values that contains both values that are smaller and values 
that arc larger than ^. The same property occurs for the representation of i. When 
considering the sum of values in the interval around | and values in the interval around 
|, the results are in an interval that contains both values that are smaller and values that 
are larger than i. However, numbers of the form 0.1 •• • can only be larger than or equal 
to ^ and numbers of the form 0.0 • • • can only be smaller than or equal to i. Thus, even 
if we know the inputs with a great precision, we must indefinitely delay the decision and 
require more precision on the input before choosing the first digit of the result: we need 
to know the inputs with infinite precision before deciding the first digit of the output. 

We solve this problem by adding a third digit in the notation. This digit provides a 
way to express that the interval given by a prefix has ^ in its interior. This new digit 
adds more redundancy in the representation. We now have three digits, even though we 
still work in base 2. We choose to name the three digits L (for left), R (for Right), and C 
(for center). 

— The digit L is used like the digit 0. If x is an infinite sequence of digits representing 
the number v, the sequence Lx represents v/2. 

— The digit R is used like the digit 1. If a; is an infinite sequence of digits representing 
the number v, the sequence Kx represents v/2 + 1/2. 

— The digit C is used with the following meaning: if x is an infinite sequence of digits 
representing the number v, the sequence Cx represents f/2 + 1/4. 

The fact that L is like 0, R is like 1, etc can also be expressed using a fimction a such 
that a(L) = 0, a(R) = 1 and a(C) = 5. With this function we can still interpret a digit 
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sequence 0.did2 ■ ■ ■ as an an infinite sum: 



^ 2* 

i=i 

From now on, we consider only numbers in the interval [0, 1] and we drop the first 
characters "0." when writing a sequence of digits. We use the same notation for a sequence 
of digits and the real number it represents. In the same spirit, we use the same notation 
for a digit d and the function it represents, the function that maps x to Last, 
we associate the digits L, R, C to the intervals [0, 1], and [i, |], respectively. These 

interval are called basic intervals. 

The redundancy of the new digit gives a very simple property: a number that can be 
written CLx can also be written LRx and a number that can be written CRx can also be 
written RLx. This property is used several times in this paper. 



3.1. Formal details of infinite sequences and real numbers 

In our formalization, we benefit from theories that state the main properties of natural 
numbers (type nat), integers (type Z), and real numbers (the type is usually written 
R, but in this article we shall write it as Rdef initions . R to avoid ambiguity with the 
"digit" R). The two integer types come with addition, subtraction, and multiplication, 
while the type of real numbers is also equipped with division. The integer types are actu- 
ally described as inductive types and the basic operations are implemented as recursive 
functions. For the real numbers, the existence of the type, two constants and 1, the 
operations, comparison predicates, and the properties of these operations (associativity, 
distributivity, inverse, etc.) are assumed. Among the assumed features, there is an axiom 
that expresses completeness, which states that every bound and non-empty subset of K. 
has a least upper bound in R. This means that whenever we exhibit a property E and 
prove that it is bound, we can construct a function that returns its least upper bound. 
This completeness axiom is inherently non-constructive. To be more precise, our work 
describes a collection of algorithms on a representation of real numbers, in some sense 
a constructive of real numbers, but many justifications of correctness, which rely on the 
axiomatized real numbers, are non-constructive. 

The axiomatization of real numbers also provides a few decision procedures. The deci- 
sion procedure field | |Delahaye and Mayero, 2001| ) solves equalities between rational ex- 
pressions, occasionally leaving proof obligations to make sure denominators are non-zero. 
The decision procedure f ourier determines when a collection of in-equations concerning 
affine formulas with rational coefficients is satisfiable. 

The type of digits is described as an enumerated type: 

Inductive idigit : Set := L I R I C. 

We provide both a numeric interpretation (named alpha) and a functional interpretation 
(named lift) to these digits. These can be defined in COQ with the following text: 

Definition alpha (d:idigit): Rdef initions .R := 
match d with L => I C => 1/2 I R => 1 end. 
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Definition lift (d: idigit) (x :Rdef initions . R) := (x+ alpha d)/2. 

The type of infinite sequences of digits is based on a polymorphic type of streams, which 
is defined as a co-inductive type using a declaration with the following form: 

Colnductive stresun (A:Set): Set := Cons: A -> stream A -> stream A. 

This definition defines stream A to be a data-type for any type A. It also defines the 
constructor Cons, with the type given in the definition. This definition is similar to a 
recursive data-type definition in a conventional functional programming language. In our 
mathematical notation, we will simply write ds instead of Cons d s. In COQ excerpts, 
we will also use the notation d::s. 

In proof systems, recursion is seldom unrestricted. In the COQ system, it is mostly 
provided as a companion to inductive and co-inductive data-structures. For inductive 
structures, the form of recursion that is provided is called structural recursion, and it 
basically imposes that a recursive function takes an element of inductive type as argu- 
ment and a recursive call can only be performed if the argument is a sub-term of the 
initial argument. For co- inductive structures, the form of recursion that is provided is 
called guarded co-recursion, and it basically imposes that a co-recursive value must be 
an element of a co-inductive type and that co-recursive calls can only be used to produce 
sub-terms of the output. More general forms of recursion are also provided, for instance 
well-founded recursion, where recursive calls are allowed only if the argument of the recur- 
sive call is a predecessor of the initial argument with respect to a relation that is known to 
be well-founded (which intuitively means that this function contains no infinite chain of 
predecessors). Well-founded recursion can actually be shown to be a special case of struc- 
tural recursion dNordstrom, 1988| [PaiTlin-Mohring.TMsl [BCTtot and Casteran, 2004| ). 

Co-recursive values need not be functions. For instance, and 1 are represented by the 
infinite sequences LLL . . . and RRR . . . , these are defined as co- recursive values with the 
following definition: 

CoFixpoint zero: stream idigit := L::zero. 
CoFixpoint one: stream idigit := R::one. 

To relate infinite streams of digits with real numbers we define a co-inductive relation. 

Colnductive represents: stream idigit -> Rdef initions . R -> Prop:= 
repr: forall d s r, 

represents s r -> <= r <= 1 -> represents (d::s) (lift d r) . 

This definition introduces both the two-place predicate represents, and a construc- 
tor, named repr, which can be used as a theorem to prove instances of this predicate. 
The statement of this theorem can be read as "for every s and r, if the proposition 
represents s r holds and the proposition <= r <= 1 hold, then the proposition 
represents (d::s) (lift d r) holds" . This relation really states that infinite streams 
are only used to represent numbers between and 1 and it confirms the correspondence 
between the digits and their function interpretations. 

An alternative approach to relating sequences of digits and real numbers is to build a 



Yves Bertot 



8 



function that maps an infinite sequence to a real value. Every prefix of an infinite sequence 
corresponds to an interval that contains all the values that could have the same prefix. 
As the prefix grows, the new intervals are included in each other, and the size is divided 
by 2 at each step, while the bounds remain rational. We defined a function bounds to 
compute the interval corresponding to the prefix of a given length for a given sequence. 
This function takes as arguments a digit sequence and a number n and it computes the 
bounds of an interval that contains all the real numbers whose representation shares the 
same prefix of length n. This function is primitive recursive in n (in other words, it is 
structural recursive with respect to the conventional representation of natural numbers 
as an inductive type). 

bounds(. . . , 0) = [0, 1] 
bounds(rfs, n+ 1) = [lift d a, lift d b] where bounds(s,n) = [a,b] 

In practice, we do not manipulate real numbers in this function, but only integers. The 
result of the function is a structure ((a, 6), k) such that the interval is 

We then define a function that maps a stream of digits to a sequence of real numbers, 
which are the lower bounds of the intervals. This function is called si_un and is defined 
by the following text, where IZR is the function that injects integers in the type of real 
numbers. 

Definition si_un (sistreeun idigit) (n:nat): Rdef initions .R := 
let (p,k) := bounds n s in let (a,b) := p in IZR a/IZR(2"k) . 

We prove that for every d, lift d is monotonic, so that si_un s is a growing sequence 
bounded by 1. All this leads us to a proof that si_un s has a limit and that this limit 
is in [0,1]. This makes it possible to define the function real_value that associates an 
infinite sequence of digits to the limit. 

We then prove that adding a digit in front of a sequence is the same as using this digit 
as a function. 

Theorem real_value_lif t : 

forall d s, real_value (d::s) = lift d (real_value s) . 

This makes it possible to show that real_value and represents follow the same struc- 
ture and to obtain the following theorem. 

Theorem represents_real_value: forall s, represents s (real_value s) . 

To complete the correspondence between the two notions we need to express that the 
relation represents is actually a function. Wc do this by expressing that the distance 
between two possible values represented by a sequence is smaller than ^ : This is proved 
by induction over n. 

Theorem represent_dif f _2pow_n : 

forall n X rl r2, represents x rl -> represents x r2 -> 
-l/(2-n) <= rl - r2 <= l/(2-n). 

It is then easy to conclude with the following theorem. 

Theorem represents_equal : forall s r, represents s r -> real_value s = r. 
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We thus have two ways to express that a given sequence represents a real number. The 
function real_value is more pleasant to use in theorem statements, but the co-inductive 
property makes proofs more elegant. Actually, all proofs of func;tion correctness presented 
in this paper are performed using a proof by co-induction based on represents, even 
though theorem statements are sometimes expressed using real_value. 

This raises a third question: given an arbitrary real number, is there a digit stream 
that represents it? The answer to this question is related to question of constructivism in 
mathematics. If you have a constructive description of your real number, then this may 
for instance be tantamount to a boolean predicate on rational numbers corresponding 
to this number viewed as a Dedekind cut (a boolean predicate that is false for every 
rational number smaller than represented real number and true for any rational number 
that is larger). We can produce the co-recursive value corresponding to any boolean 
predicate using a co-recursive function. Here is an example, where the rational numbers 
are viewed as pairs of integers (p (a, b) = true means that the real number of interest 
is smaller than or equal to |): 

CoFixpoint stream_of _cut (p:Z*Z->bool) : stream idigit := 
match p (1, 2) with 

true => R: : strecim_of _cut (fun r => let (a,b) := r in p (a+b, 2*b)) 
I false => L: :stream_of _cut ifxm r => let (a,b) := r in p (a, 2*b)) 
end. 

Alternatively, the real number of interest can be given by a Cauchy sequence of rational 
numbers and a constructive proof that it satisfies the Cauchy criterion. One way to 
describe Cauchy sequences is to fix a function h from nat to R, with its limit in when 
the argument goes to infinity. A Cauchy sequence may then be given by a function / 
from Z to Q and the Cauchy criterion may be given by a monotonic function g from Z 
to Z such that 

Vn m p,g{n) < m A g{n) < p |/(m) - f{p)\ < h{n) 

To construct the infinite list of digits for a given Cauchy sequence, we simply need to 
repeat the following process: 

1 compute the n first elements of the stream, this actually gives an interval of width 
^ , compute the lower bound b of this interval, 

2 find the least n such that h{n) < 2?n^, 

3 compute f{g{n)), we know that the distance between this value and the sequence's 
limit is less than ^^_^3 , 

4 compute the value a = 2'"-{f{g{n)) — 6), this value is in [0,1] by an invariant of the 
recursive process, 

5 if a < |, then we know that for every m> g{n), f{m) < | and we choose the n-|- 1th 
digit to be L, 

6 if I < a < |, then we know that for every m > g{n), j < f{m) < |, and we choose 
the n + 1th digit to be C, 

7 if I < a, then we know that for every m > g{n), ^ < f{m), and we choose the n+lth 
digit to be R 
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We contend that this technique gives a constructive process to associate a sequence of 
digits to a sequence of rational numbers associated with a constructive proof that this 
sequence satisfies the Cauchy criterion. 

When implementing this recursive process as a co-inductive function, we propose to 
represent rational numbers with pairs of integers and to replace comparisons of rational 
numbers with comparisons of integers. It is also more convenient to restrict ourselves to 
the case where h{n) = In the recursive process, wc do not need to keep the list of the 
n first digits, we only need to know n and the lower bound of the represented interval, 
these are given as arguments to the co-recursive function: 

CoFixpoint streain_of _cauchy (f : Z -> Z*Z)(g: Z -> Z) 
(n:Z)(b:Z*Z) : stream idigit := 
let (vn,vd) := f(g n) in 
let (bn.bd) := b in 
let (d, r) := 

if is_smaller (8*2"n* (vn*bd-vd*bn) ) (3*vd*bd) then 
(L,b) 

else if is_smaller (8*2"n*(vn*bd-vd*bn)) (5*vd*bd) then 

(C, (4*2~n*bn+bd,4*2"n*bd)) 
else (R, (2*2"n*bn+bd,2*2~n*bd)) in 
let (new_bn, new_bd) := r in 

d: : stream_of _cauchy f g (n+1) (new_bn, new_bd) . 

The functions stream_of _cut and stream_of _cauchy are only given here to show the 

feasibility of connecting streams of digits with the usual notions of real number con- 
structions, but more efficient ways to handle rational numbers should be used if these 
functions were to be used effectively, for instance least common denominators should be 
computed between fraction numerators and denominators. 

For an arbitrary real number between and 1 given non constructively, but known 
by its binary representation (an infinite sequence of and 1 digits), this real number is 
simply represented by the corresponding infinite stream where is replaced with L and 1 
is replaced with R; however, the fact that the real number is given non constructively is 
reflected by the fact that we can't write a co-recursive function that produces the stream. 

3.2. Addition 

It is well-known''' that adding two infinite sequences of redundant digits can be described 
as a simple automaton that reads digits from both inputs and produces digits at every 
recursive call. Two approaches can be taken: either this automaton is understood as a 

program that keeps a carry as it processes the inputs, or it can be viewed as a program 
that performs a little look-ahead before outputting the result and processing the rest, 

t P. Martin-L6f suggested to the author that Cauchy had aheady devised an algorithm for addition in 
a similar representation. Di Giaiiaiitonio refers to Cauchy and Leslie, but the reference to Cauchy's 
work is wrong and the reference to Leslie could not found at the time of writing. 
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maybe with a slight modification of the first digit in both inputs. The algorithm we 
describe follows the second approach. 

Our algorithm has two parts. The first part is a function that computes the arithmetic 
mean of two values (in other words, the half-sum). The second part is a function that 
computes the double of a value. The first algorithm maps two real numbers in [0, 1] to a 
value in [0, 1]. The doubling function only returns a meaningful value when the input is 
smaller than or equal to i. 

For the half-sum, the structure of the algorithm is as follows: if the inputs have the 
form did2X and d^d^y, then the algorithm outputs a digit d and calls itself recursively 
with the new arguments d^x and dey. Written as an equation, this yields the following 
formula: 

half _sum((ii(i2.T, fi3(i4?/) = (i(half _sum((i5a;, rfgy)). 

Here is an example, suppose that di =L and = R, in this case we can choose d = C 
and ^5 = d2 and de = d4, because the following equalities hold, using the interpretations 
of digits as functions: 



half _sum(Ld2a;, Bd4y) 



Ld2X + Rd4y 



d-i'x I rf.ij; I _L 
2 ^ 2 ^ 2 

2 

d2X d^y 1 

2 I 



2 4 
= C(half _sum((i2a;, (i42/)). 

In this case, it is not necessary to scrutinize d2 and d^ to decide the value of d and the 
arguments for the recursive call. The equation can be re-written as 

half _sum(La;, Ry) = C(half _sum(x, y)). 

Here is a second example where d^ and d^ are modified with respect to d2 and d^. We 
suppose d\ = C, d2 = L, and d^ = L. In this case, the following equalities hold: 



half _suin(CLa;, Ldiy) 



1 I djJl 

4^2 



2 

I 1 I (Uy 
2 2 



2 

= L(half _sum(Ra;, ^4?/)) 

In this case, it is not necessary to scrutinize ^4, but the value d^ is modified with respect 
to rfe- The equation can be re- written as 

half _sum(CLa;, Lj/) = L(half _suiii(Ra;, t/)). 

If we had designed the algorithm to scrutinize two digits in each input, there would be 
81 cases, but since some cases can be handled with a scrutiny of only the first digit in 
each input, or only one digit in one of the inputs, the number of cases is brought down 
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to 25 cases. Moreover, the exact behavior of the algorithm in each case can be found 
automatically with the help of proof search procedures. 

3.3. Automatic generation of function code 

We use the proof engine to actually construct the half-sum function by making this 
program use its proof search facility to construct a term with the right type, which 
should be 

stream idigit -> stream idigit -> stream idigit 

We have enough control on the proof search mechanism to express that the value of 
this type that wc want to construct should be a co-rccursivc function and that it should 
analyze the first digit of the two input streams. This is simply done by stating that a 
case analysis on these arguments should be performed. Doing a case analysis on the first 
digit of the first input yields three arguments, doing a case analysis on the first digit 
of the second argument also gives three cases, so that there arc at least nine cases to 
consider. Some cases are easily solved directly, by simply finding an output digit that 
makes the addition correct. For instance, if the two inputs are dx and dy (in other words, 
they have the same initial digit) then the result should be rf(half _siim x y). This because 
the formula holds: 



When no output digit can be computed to make the formula work directly, more in- 
formation is gathered from the inputs by imposing more case analysis. When this case 
analysis is performed, we look at possible values of the second digit of one of the inputs 
and we decide if we have enough information to decide what the first output digit should 
be. This decision is taken by performing some interval reasoning. If two digits of one 
of the inputs are fixed, then this input belongs in an interval of length \, adding this 
interval with an input for which only one digit is known, this gives an interval of length 
|, taking the half of this gives an interval of length |. If the lower bound of this interval 
is 0, , i, I;, ... then wc know what the output first digit can be, but if the lower boimd 
of this interval is ^ (this happens when one of the inputs is LCx and the other is Cy) , 
then the upper bound is ^ and we cannot conclude because this interval may contains 
values lower than j (which should not start with a C or a R) and values larger than i 
(which should not start with a L): for these cases an extra case analysis on the second 
input is required. 

When the first digit of the output is fixed, we still need to decide what the first digit 
of the arguments to recursive calls will be. This may include a change with the initial 
second digit of the input. This difference is often related to the equivalence between LR 
and CL prefixes on the one hand and between RL and CR on the other hand. 

For instance, if the first input has the form CLa; and the second input has the form 
Li/, the function can return L(half _sum(R.x y)), because the half sum of the two inputs 
is equivalent to the half sum of LRx and Ly. 



2 
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To determine the first digits of inputs in recursive calls, we must first respect an 
important rules: variables that represent sub-streams should appear behind the same 
number of digits in the input pattern and in the output pattern. This rules comes from 
the fact that the real number represented by these variables is divided by 2 every time a 
digit is added in front of it. If we want the half-sum equation to be satisfied we must make 
sure that the number of divisions by 2 is preserved between the inputs and the output. 
Thus, we can only prescribe the first digit of one the arguments to the co-recursive call 
of half _suin if the corresponding input had two digits in the input pattern. 

To determine what first digit can be used for the recursive call on an argument, the 
proof search procedure compares the lower bound of the output as prescribed by its 
first digit to the lower bound of possible results that we can predict from the shape of 
the input patterns. In the example, the lower bound of the output as prescribed by the 
first digit L is 0, the lower bound predicted from the half-sum of CLx and Ly is |. The 
discrepancy must be resolved by making sure that the sum of all the digits appearing as 
head of recursive call arguments adds up to i (which does fit with a target i since we 
are computing a half-sum and place the output's first digit L in front). Here there is only 
one digit available, and we can only choose its value to be R. 

Although the half-sum function is obtained by mechanical means to be correct by 
construction, its type is only 

stream idigit -> stream idigit -> stream idigit 

This type does not express what the function does. We need to add a theorem to state 
that it has the right behavior with respect to the real numbers represented by the inputs 
and outputs. The statement has the following shape. 

Theorem half _sum_correct : 

forall X y u V, represents x u -> represents y v -> 
represents (half_sum x y) ((u + v)/2). 

The proof of this statement relies on a proof by co-induction: we assume that the state- 
ment is already satisfied for any output stream that is a strict sub-stream of the current 
output and we show that this is enough to prove the result for the current output. The 
proof analyzes the behavior of the half _sum function and explores all the 25 cases that 
were found at the time the function was constructed. In all cases, it is a simple matter 
to verify the equality between the algebraic formula corresponding to the half-sum of the 
inputs and the output pattern present in the half _sum function, using a tactic named 
field ( |Delahaye and Mayero, 200 1| ) that solves equalities between rational expressions 
in a field; a second statement that needs to be verified is that the half-sum of the inputs 
does belong to [0,1] if the two inputs do, this is easily done using a tactic named f ourier 
that solves affine comparisons between real values with rational coefficients. 

We believe that our definition technique can easily be reproduced for different sets of 
digits or for other simple binary operations like subtraction. 
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3.4. Multiplication by 2 

We also need to provide a function to multiply the output of half _sum by 2. This function 
is based on the following remarks. 

— The double of a number of the form Lx is simply x, 

— The double of a number of the form Ra; is either 1 or outside the interval [0,1], 

— The double of a number of the form Cx is a number of the form Kx' where x' is the 
double of X (hence the algorithm exhibits a co-recursive call). 

Here is the formal definition: 

Cofixpoint mult2 (v:stream idigit) : stream idigit := 

match V with L: :x => x | C: :x => R: :inult2 x | R: :x => one end. 

The correctness of this function is expressed with the following statement: 

Theorem mult2_correct : forall x u, 

<= u <= 1/2 -> represents x u -> represents (mult2 x) (2*u) . 

Please note that this theorem explicitly states that the result value is specified correctly 
only when the input is smaller than ^ . 

3.5. Subtraction 

In this section wc discuss several approaches to subtraction. One first approach uses a few 
intermediary functions. The first intermediary function mimics the opposite function. Of 
course the opposite function cannot be defined from [0,1] to [0,1], but we can define the 
function that maps a; to 1 — x. Here is the general definition, where we name the function 
ininus_aux: 

minus_aux(L(a;)) = R(ininus_aux(a;)) 
minus_aux(C(x)) = C(minus_aux(x)) 
minus_aux(R(a:)) = L(minus_aux(x)) 

These equations are justified through simple computations. For instance, the last equation 
is justified with the following reasoning steps: 

X 1 

minus_aux(R(a;)) = ~ ^2 2"^ 

1 X 

2 ~ 2 

= i(l-x) 
2^ ' 

= L(minus_aux(.T)) 

Combining ininus_aux with addition, we can easily compute the binary function that 
maps X and y to 1 + x — y. Of course, this function returns a meaningful result only when 
X is smaller than y. 

Alternatively, we can combine minus_aux with half _sum to have a function that maps 
X and y to , Now, if we really want to have a subtraction, we can remove the ^ 
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offset. We use another auxiliary function, which we name minusJialf and is defined by 
the fohowing equations: 

minusJialf (Rx) = La; 
minusJialf (Lx) = zero 
minusJialf (Cx) — L(minusJialf (x)) 

The first of these equations is trivial to justify. The second is justified by the fact that 
the only value inside [0, ^] for which x — ^ belongs to [0, 1] is ^ and the result is in this 
case. The third equation is justified by the following reasoning steps: 

X I I 
minusJialf (Cx) = o' + T ^ 



X 1 
2 ~ 4 




= L(minusjialf (x)) 



4. Parameterized affine operations 

In this section, we study another approach to addition, with the encoding of a more 
general function that computes affine formulas in two real values with rational coefficients. 
More precisely, we want to compute the value 

a b c 
a' b' d 

When a, 6, c are non-negative integers and a', c' are positive integers (a, 6, c may be 
null, but not the others), and x and y are real numbers, given as infinite sequences of 
digits. 

The interpretation of digits as affine functions (using our function lift) makes them a 
special case of what Edalat and Potts call Linear Fractional Transformations ( |Edalat and Potts, 1998| ). 
They actually show that a more general form of two argument transform can be pro- 
gramed on this form of real number transformation, namely the computation of the 
following transform, called a Mobius transform, where a, 6, etc. are integers: 

axy + bx + cy + d 
exy + fx + gy + h' 

Studying only affine formulas correspond to restricting the general study proposed by 
Edalat and Potts to the case of where e, /, and g are 0. A good reason to study separately 
the restricted case is that the formal proofs stay inside the realm of proofs about equali- 
ties and comparisons of affine formulas with rational coefficients, a realm for which auto- 
matic proof tools exist at the time of writing this article, while verifying the correctness 
of the general Mobius transform requires incursions into the realm of proofs about equal- 
ities and comparisons of polynomial formulas, a domain for which proof procedures are 
only under development ( |Mahboubi and Pottier, 2002) [McLaughlin and Harrison, 2005| 
IMahboubi, 2006| ). 
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4.1. Main structure of the algorithm 

Choosing the digits of the result is based on the following remarks: 

1 Even without observing a; and y, we already know that they are between and 1. 
The result lies between the extrema 

c ab'c' + a' be' + abc' 
d ■ 

2 if the extrema belong to the same basic interval, it is possible to produce a digit and 
perform a co-recursive call with a new afhne formula. In this sense, the output does 
not depend on reading more digits from the input and the algorithm can be described 
as a streaming algorithm in the sense of ( [Gibbons, 2004| ). 

3 If the extrema are badly placed, we cannot choose an interval associated to a digit 
that is sure to contain the result. In this case, we scrutinize x and y and observe their 
first digit. As a result, we obtain a new estimate of the interval that may contain the 
result, whose size is half the previous size. We can then perform a recursive call with 
a new affine formula. In the long run, we are forced to get to a situation where the 
extrema are inside a basic interval and a co-recursive call can be performed. In fact, 
this condition is guaranteed as soon as the distance between extrema is shorter than 
1/4. 

Let us study two examples. In the first example, suppose that the property ^ > 1/2 
holds. We know that the result is larger than 1/2 and we can produce a R digit. The 
following computation takes place: 

a b c /„ , a b c , , 

-^+T;y+- = H2x{-x + -y +-)-!) 
a' b' d a' b' d 

,2a 2b 2e-d, 
a' b' d 

There is a recursive call with a new affine formula, where all the coefficient are positive 
integers or non-negative integers as required. 

In a second example, suppose that the properties x — Lx' and y — Ry' hold. The 
following computation can take place: 

a b c , ^\ ,by' + l c 

_ a , b , bd + 2b'c 

~ 2^^ ^ 2&^^ ^ 2b'd 
Here again, we can have a recursive call with a new affine formula, no digit has been 
produced (therefore the recursive call cannot be a co-recursive call) but the distance 
between the extrema in the new formula is a /2a' + 6/26', the half oi a/a' + b/b', which 
was the distance between extrema for the initial affine formula. 



4.2. Formal details for affine formulas 



When providing the formal description of the recursive algorithm for the computation of 
affine formulas, we need to pay attention to the following two important aspects: 
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1 The function needs to be partial, because we must ensure the sign conditions on the 
coefficients of the affine formula, 

2 Not all recursive calls are acceptable co-recursive calls, because recursive calls after 
the consumption of digits from the two input streams are associated to no production 
of an output digit. 

We define a record type named af f ine_data that collects the eight elements of an affine 
formula and a predicate named positive_coef f icients that states that the coefficients 
satisfy the sign conditions. 

The computation is then represented by a main function called axbyc to compute the 
affine formula. This function has a dependent type: it takes as first argument an affine 
formula and as second argument a proof that its coefficients satisfy the sign conditions. 

axbyc: forall x: af f ine_data, positive_coeff icients x -> stream idigit . 

This function is defined as a co-recursive function. The constraints on recursive program- 
ming impose that this function can only perform the recursive calls that are associated 
to the production of a digit in the output (phase El in the previous section). We need 
to define an auxiliary function, not a co-recursive one, which takes charge of the recur- 
sive calls that are only associated to the consumption of digits from the input streams 
(phase El in the previous section). 

The auxiliary function is named axbycjrec. It takes as arguments an affine formula and 
a proof that it satisfies the predicate positive_coeff icients. It returns an equivalent 
affine formula, one for which the decision of producing the next output digit can be 
taken. This is represented by the fact that output of this function is in a type with 
three constructors, called caseR, caseL, or caseC. Each constructor contains as first 
field the new affine formula, as second field a proof that this new formula satisfies the 
sign conditions. The next field (for the constructors caseR and caseL) or the next two 
fields (for the constructor caseC) express that the right interval conditions are satisfied 
to output a digit, The last field is a proof that the new affine formula is equivalent to 
the initial one. 

The recursive structure of the function axbycjrec is based on well-founded recursion. 
More precisely, it relies on the fact that the distance between extrema decreases at 
each recursive step. This can be translated into an integer formula that decreases while 
remaining positive. When this integer formula is 0, we can prove that one of the interval 
conditions to output a digit is necessarily satisfied. 

Two other collections of auxiliary function perform the elementary operations. Func- 
tions named prod_R, prod_L, and prod_C perform the coefficients transformations that 
should be performed after producing an output digit. For instance prod_R maps the affine 
formula 

a b c 
a' b' c 

to the formula 

2a 2b 2c - c' 

— a; + Try + — — 

a' b' c 
as we already justified in the first example of the previous section. 
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The lemmas named A.prod_R_pos, A.prod_L_pos, and A.prod_C_pos ensure that the 
functions prod_R, prod_L and prod_C preserve the sign conditions on the coefficients, 
respectively. These lemmas rely on the interval conditions produced by axbyc_rec. For 
instance, for A.prod_R_pos the extra interval condition is c' < 2c. In our description of 
the co-recursive function, this information is transferred from axbyc jrec to A . prod_R_pos 
with the help of a variable named He. 

With these auxiliary functions, the main function can be given a simple structure: 

CoFixpoint axbyc (x : af f ine_data) 

(h:positive_coef f icients x):streain idigit := 
match axbyc_rec x h with 
caseR y Hpos He _ => 
R: : (axbyc (prod_R y) (A.prod_R_pos y Hpos He)) 
I easeL y Hpos _ _ => 

L:: (axbyc (prod_L y) (A.prod_L_pos y Hpos)) 
I easeC y Hpos HI H2 _ => 

C:: (axbyc (prod_C y) (A.prod_C_pos y Hpos H2)) 
end. 

With the help of the function real_value we can also define a function af _real_value 
that maps any affinc formula represented by an element of aff ine_data to the real 
number it represents. This function is used to express the correctness of our algorithm 
to compute the affine formula: 

axbye_eorreet : 

forall X, forall H :positive_coeff icients x, 
<= af _real_value x <= 1 -> 
real_value (axbyc x H) = af _real_value x. 

This proof is based on a lemma that is proved by co-induction: 

axbye_eorreet_aux : 

forall X : af f ine_data, forall H : positive_eoeff icients x, 
<= (af _real_value x) <= 1 -> 
represents (axbyc x H) (af _real_value x) . 

This algorithm is interesting because it provides ways to add two real numbers, to multi- 
ply them by rational numbers, and to encode rational numbers as real numbers. Having 
formalized this algorithm may make the direct implementation of addition described ear- 
lier seem useless. It is not useless, because the direct implementation of addition makes 
no use of dependent types, proof arguments, or well-founded recursion. As a result, the 
direct addition can be executed directly inside the theorem prover using its internal re- 
duction mechanism, while the affine formula computation can only be executed outside 
the theorem prover as extracted code. Algorithms that can be reduced inside the proof 
system may play a role in reflection-based proof procedures ( [Boutin, 1997| ). 



Affine functions and series with co-inductive real numbers 



19 



5. Computing series 

A series is an infinite sum of values. Knowing how to compute series can help in computing 
famous constants like e (Euler's number) or tt and to implement the multiplication of 
two real numbers. 

5.1. Main structure oj the algorithm 

We consider the problem of computing values of the form X^^q a, where the terms 
are real numbers. 

Studying series is very close to studying converging sequences, since it is enough to 
consider the sequence u„ = X^"^o '^i- Each element of the sequence can then be computed 
as a finite combination of additions. 

Computing the first p digits of the limit means computing the limit up to a precision 
of 2~'P. If we know that a given element u„ is closer than 2"^^'+^^ to the limit and if 
we can compute u„ to a precision better than 2^^+^, then we are able to compute the 
limit to the required precision. In other words, if we know that | | is smaller than 

2-(p+i), and we know a value v whose distance to Y^^=o 1^^^ than 2~(p+^\ then we 
also know that the distance of v to is less than X^i^o '^i- fact, wc need to use slightly 
stronger precisions, because an interval of length 2~p rarely fits in one of the intervals 
representable by finite sequences of digits. Still, this approach shows that we can avoid 
considering the whole infinite sum to produce the first digit of the output. 

Wc restrict our study to series whose convergence is described constructively by a 
function ji that satisfies the following properties: 

oo 

Vm. n < m ^ I \^ ai\ < iJ,{n) lim fj,{n) = 

i=m 

We actually formalize the computation of a function / that has the following informal 
specification: 

oo 

f{x,y,n) = X + y x^tti, 

i—n 

where a; is a real number, y is a rational number, and n is a natural number. Intuitively, 
y represents the inverse of the precision that is reached in the computation {y = 2^). 
If we know that y x (i{n) < and we know x to a precision of three digits then we 
are able to choose the first digit of a; + yJ2^n ^i- ^® then perform the following 
computation 

f{x,y,n) = df{2x ~ a{d),2y,n), 

In most cases, we also have x = dx' for some x', and the value 2x — a{d) is simply 
represented by x' . If y x /x(n) > we compute a new value (j){y,n) such that y x 
IJ,{4'{y,n)) < jq. This value is bound to exist because /i converges to at infinity. We 
can then proceed with the following step. 

<^(?/,n)-i 

f{x,y,n) = f{x + y X ^ ai,y,(f){y,n)) 

i=n 
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In practice we first compute (t){y,n), then we compute v = x + y X]fi^'"^ ^ o,i by 
repeated binary additions. Let us assume that p is defined as y x X]i^0(yra)^«- ^'^^ 
number we want to compute is u + p and we know that \p\ < We then perform the 
following case analysis: 

1 If u has the form Rw', and v' has the form R?;", C7/', LC?;", or LRw", wo can deduce 
that V > 9/16, therefore v + p > ^ and the first digit of the result can be R. The 
result is K{f{v',2y,(j){y,n))). 

2 li V has the form Cw', where v' has the form Cu", LCw", LR?;", RL?;", or RCw", then we 
are certain that jq <v < therefore \ < v+p < |, the result is C{f{v', 2y, (j){y, n))). 

3 If I) has the form Lv', where v' has the form Lv", Cv", BLv", RCv", then we are certain 
that V < jq and v + p < ^ the result is L{f{v', 2y, (f){y, n))). 

4 a V has the form RLLv" , then ?; could also be represented using CRLf" and this case is 
already considered above. The same goes for the cases LRR, CLL, and CRR, using CLR, 
LRL, and RLR as respective alternatives. 

The number (j){y,n) is chosen so that y x p{(f){y,n)) < jq because is the shortest 
distance between the boimds of the intervals for CLC and C, RLC and R, or LRC and 
L. Computing (j){y,n) and X^fi^'"''"^ ai depends on the series being studied. Because 
recursive calls to / always have 2y and </>(?/, n) as arguments, we can also assume that 
the invariant y x /x(n) < | is maintained through recursive calls. 

5.2. Series with positive terms 

When we know that the a, terms are all positive, we do not need to use ^ to bound the 

infinite remainder of the series. The computation technique can be simplified. 

We first compute (/>(y,n) so that y x p{<j)[y,n)) < | and v = x + X^fl^'"^""^. In what 
follows, let p be defined as y x J2il,p{y n) know \p\ < |. We perform the following 

case analysis: 

1 if V has the form Kv', we are sure that the result is larger than or equal to 1/2, the 
result is R(/(u', 2y, <j){y, n))). 

2 if i> has the form Cv' but not CRw', we can deduce v G [3, |] and v + p £ [3; |]> the 
result is C{f{v',2y,(p{y,n))). 

3 ii V has the form Lu', but not LRv', we can deduce v € [0, |] and v + p G [0, ^]. The 
result is L{f{v',2y,(j){y,n))). 

4 ii v has the form CRv" or LRv", we can use the equivalences with RLv" and CLv", 
respectively, to switch to one of the previous cases. 

For positive series, there is also an invariant: y x /^(n) is always smaller than j. This 
invariant plays a role in the correctness proofs. 

5.3. Formal details of positive series 

We programmed the general treatment of positive series in a higher-order function that 
is easy to re-use from one series to the other. 
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Definition series_body (A:Set) 

(f : stream idigit -> A -> stream idigit) 
(x: stream idigit) (a:A): stream idigit := 
let (d,x') := V in 
match, d with 

R => R: :f v' a 

I L => match v' with R::v" =>C::f (L::v") a I _=>L::fv' a end 
I C => match v' with R: :v" => R: :f (L::v") a I _ => C::f v' a end 
end. 

The type A should make it possible to compute the values y and n used in our informal 
presentation. This will be made clearer in the next examples. 



5.4. Application to computing e 

The number e is defined by the formula 



oo ^ 
k=0 



Of course, this number is larger than 2, but we only want to compute its fractional part, 
so that we actually compute Xlfc^2 IT- following properties are easy to obtain, by 
remembering that n\n^~'^ < k\ for every k >n, therefore 



o< 



£-fc! (n-l)!(n-l)- 



k=n 

For this series, we choose /u(n) to be the value ^n-iy.{n-i) ■ soon as 2 < n, we have 
/i(n + 1) < ^i^. This implies the following property: 

Vn, y,0 < y A2 < n Ay X fi{n) < ^ ^ n) <n+l. 

Thus, it is never necessary to absorb more than one term from the infinite sum into x at 
each co-recursive call. 

the type A that appears in our use of series_body is a triple type. The triple given 
as argument has the form {y, n, 6), where is the precomputed value 9 = {n — 1)!. This 
invariant is maintained through recursive calls so that factorials are not recomputed from 
scratch each time. Here the function is given by the formula 

^^"^ " (n-l)!(n-l) " e X (n- 1)' 

and computing (/>(y, n) is easy, because we know that this value is always n or n + 1. To 
decide whether n) = n, we simply need to compare y x with i, in other words 

to compare 8y with /x' = i = 6[n — 1). In the following code, rat_to_stream builds the 
infinite sequence of digits for a rational number given by its numerator and denominator. 

CoFixpoint e_series (x: stream idigit) (s :Z*Z*Z) : stream idigit := 

let (aux, theta) := s in let (y, n) := aux in let mu' := theta*(n-l) in 
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let (v, phi, theta') := 
if Z_le_gt_dec (8*y) mu' then 
mk_triple x n theta 

else 

let theta' := mu'+theta in 

mk_triple (x+(rat_to_stream y theta' )) (n+1) theta' in 
series_body _ e_series v (2*y, phi, theta'). 

To express that this function computes correctly the scries, we rely on a predicate 
inf inite_sum which takes a function / from Z to R and a value in R as arguments 
and means X^^q /(*) converges and the limit is v. The correctness statement has the 
following shape: 

Theorem e_correctl : 
forall V vr r y n, 

4 * y <= fact(n-l) * (n-1) -> 2 <= n -> 1 <= y -> 
represents v vr -> 

inf inite_suin (fun i => l/IZR(fact (i+n))) r -> 
vr + (IZR y)*r <= 1 -> 

represents (e_series v (y, n, fact(n-l))) (vr+(IZR y) *r) . 

This statement is obfuscated by the simultaneous use of two types of numbers and the 
function IZR is the natural injection of integers into the type of real numbers. In this 
statement the formula fact (n-1) * (n-1) represents the inverse of /z(n) and the formula 
4 * y <= fact (n-1) * (n-1) corresponds to the invariant of the series. Note that this 
theorem expresses that the scries is correctly computed only if the scries really converges 
to a value that is smaller than or equal to 1. The proof that the series converges has to 
be done independently. 

We initialize the recursive computation with x = ^, y = 1 and n = 3, so that the 
invariant on y x /i(n) is satisfied. 

Definition number_e_minus2 : streajn idigit := 

e_series (rat_to_stream 1 2+rat_to_streain 16) (1, 4, 6). 

The correctness theorem then makes it possible to obtain the following statement; 

Theorem e_correct : 

infinite_sum (fun i => l/IZR(f act (i+2) ) ) (real_value nijmber_e_minus2) . 

This statement really means IT = number_e_minus2. 

The value number_e_minus2 can then be used to construct rational approximations of 
e — 2, using the bounds function. Actually, given n we compute (a, b, k) so that 



For n = 320, this computation takes approximately a minute with the standard version 
of COQ-'- . These functions can be used in proof procedures to compute approximations 

* Coq version 8.0pl2, Intel Pentium(R) M ITOOMhz. 
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of e. The code can also be extracted both to Ocaml and to Haskell. Running the Ocaml 
extracted code, we can compute 2000 redundant digits of e — 2 in about a minute. 

5.5. Multiplication as a special case of series 

When u is the infinite sequence did2 ■ ■ ■ and then uu' is a series: 



This is a series where all terms are positive. Moreover, two simplifications can be made 
with respect to the general approach. First, while y is multiplied by 2 at every recursive 
call, Oi contains a divisor that is also multiplied by 2, so that the two multiplications by 2 
cancel out. Second, it is reasonable to simply consume one element from the infinite sum 
at each recursive call, without scrutinizing the value of u' . If this approach is followed, 
the argument y is not necessary anymore: only the digits di and u' are needed. We can 
re-use the general function series_body in the following manner: 

Definition sum_mult_d (d:idigit) (u v: stream idigit) := 
match d with L => u I C => u+L: :L: :v I R => u+L: :v end. 

CoFixpoint mult_a (x:streajn idigit) (p : stream idigit*stream idigit) 
: stream idigit := 
let (u,v) := p in 

match u with d::w => series_body _ mult_a (sum_mult_d d x v)(w,v) end. 

The function mult_a x {u, u') computes x + uu' only when uu' < j (here again, we see 
the invariant of the general approach). To obtain multiplication for the general case, we 
divide the second operand by 4 before the multiplication and we multiply the result by 
4. Here is a naive implementation: 

Definition mult (x y:stream idigit): stream idigit := 
mult2(mult2(mult_a zero (x,L: :L: :y) ) ) . 

The following theorem can then proved and verified formally: 

mult_correct 

: forall (x y: stream idigit) (vx vy: Rdef initions . R) , 

represents x vx -> represents y vy -> represents (x*y) (vx*vy) 

In this statement the notation x * y represents our multiplication as an operation on 
infinite digit streams, while the notation vx * vy represents the multiplication of real 
numbers, as they are axiomatized in the COQ system. 

It turns out our approach of multiplication, based on series, yields an implementation of 
multiplication that is very close to the implementation in CCiaf faglione and di Gianantonio, 2000| 

We can also try this multiplication directly inside COQ, for instance, we can compute 
(e — 2)'^. With the current standard version of COQ no effort is made to exploit possible 
sharing in the lazy computation of values, so that the same value may be computed 
several times. For this reason, we cannot compute this number to a high precision as 
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easily as for e — 2. For example, our few experiments showed that it takes approximately 
10 seconds to compute an approximation with an accuracy of ^ digits and a minute to 
compute the approximation with an accuracy of ^str ■ 



6. Conclusion 

The work described in this paper is available on the Internet at the following address: 

\protect\vrule widthOpt\protect\href {http : / /hal . inria. f r/ inria-00001171}{http : //hal . inria.fr/inr: 

This is both an extension and a departure from the work of Ciaffaglione and di Gi- 
anantonio. We chose to work with an extra digit that is interpreted as ^ instead of —1. 
We believe algorithms are easier to design in this approach and easier to prove correct, 
but this approach is less well adapted to expanding to the full real line. 

We believe the choice to rely on an existing axiomatization of real numbers was in- 
strumental in making the proofs in this experiment quicker to obtain. In particular, we 
relied on a collection of automatic decision procedures for equalities between polynomials 
and sets of in-equalities between affine formulas. This axiomatization relies on classical 
mathematics, and we don't know in which respect the decision procedures rely on these 
these classical axioms and in which respect they could be re-implemented in construc- 
tive mathematics. Admittedly, we only provide algorithms to compute representations 
of real values that are constructively definable. We believe that our experiment should 
be reproducible in the context of constructive mathematics, but it concentrated on the 
presentation of real numbers that provided the most complete theory of functions. The 
question of constructive or non-constructive mathematics was secondary in this experi- 
ment. 

It is characteristic that the definition of series appears to be more basic than multipli- 
cation, but this is simply due to the fact that the digit-based representation of numbers 
already is naturally intcrpretable as a series and this structure is preserved by multipli- 
cation thanks to distributivity. 

Now that we have a multiplication for our representation of real numbers, we can 
consider the task of implementing other functions, like division and analytic functions. 
For division, we expect to use a method of range reduction: to compute ^ we should 
compute ^ or I — fc where k should be chosen so that the result is inside [0, 1] but 
finding the right value for k is only possible when we have a constructive way to prove 
that y is non zero. 

An alternative approach to division is provided if we generalize our approach to affine 
formulas to consider Mobius transforms (already studied in ( |Edalat and Potts, 1998| )): 

, . axy + bx + cy + d 

{x,y) I— > . 

exy + fx + gy + h 

However, we suspect that the proofs of correctness for these transformations are less easy 
to automatize, because they do not rely only on affine formulas (which are easily treated 
with the Fourier-Motzkin decision procedure. 
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One of the interesting features of this experiment is that some series can be computed 
directly inside the theorem prover. This is an important feature, if a proof requires that 
we produce an accurate approximation of a series value. 

Concerning the computation of mathematical constants, we already have all the in- 
gredients to compute tt using a Machin formula (we used the formula j = arctan(i) + 
arctan(i) which can easily be proved and computed as the sum of two series). 

Several questions will be studied in the future. First, the choice of a digit set is arbi- 
trary. We already experimented we a digit set that contains only two digits, corresponding 
to intervals [0, |] and [|, 1]; although some functions seem to be simpler (because there 
are less cases to consider, for instance when considering afhne formulas) other problems 
arise because equalities between patterns do not exist for short patterns (thus we cannot 
as easily mimic the equality CL — LR). di Gianantonio also studied a two-digit repre- 
sentation and he proposes to work with intervals whose length is based on the golden 
number ( |di Gianantonio, 1996| ). However, the price to pay is that we now need to solve 
a polynomial system of degree 2 to determine whether a given value belongs to one of 
the basic intervals, again formal proofs in this setting are made more complex because 
we are not in the real of affinc formulas anymore. 

In the long run, we wish to choose digit sets that are closer to the bound integers 
that are usually handled in computers, so that regular integer addition, subtraction, 
multiplication, and comparison or even bitwise operations can be used to establish the 
basic relations between various digits. 

The second important question is related to the efhciency of co-recursive computation 
inside the theorem prover. While recent evolutions of the COQ system brought drastic 
improvements in the computation of recursive functions, it is not certain that the co- 
recursive question is as well treated. In particular, lazy computation is needed to achieve 
reasonable speeds, while the current version of COQ may only implement call-by-name. 
This question is particularly difficult because the reduction mechanism also needs to 
retain the property of strong normalization for non-closed terms. 

Eventually, we hope to achieve the development of a reasonably efficient and formally 
verified library for mathematical computation. We believe this will be a stepping stone 
for more ambitious projects like the Flyspeck project plales, 2000| [Hales, 2004| ). 
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